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We derive a model that describes the motion of a Brownian particle in a system which is dominated 
by gravitational forces. An example of such a system is a massive black hole immersed in a cluster 
of stars. We compute the dispersion in the position and velocity of such a black hole, and examine 
whether it achieves a state of equipartition of kinetic energy with the stars. This problem has been 
considered before only for stellar systems with an isothermal Maxwellian distribution of velocities; 
here we study other examples and confirm our calculations with N-body simulations. In certain 
cases, depending on the stellar distribution function, the black hole can acquire a steady state 
kinetic energy which is very far from equipartition relative to the stars. 

PACS numbers: 95.10.Ce, 45.50.Pk, 05.40. Jc, 98.10.+Z, 98.62. Js 



Brownian motion in stellar systems was studied by 
Chandrasekhar who recognized that the force act- 
ing on an object in a stellar system consists of two 
independent contributions: one part, which originates 
from the "smoothed-out" average distribution of mat- 
ter in the stellar system, varies slowly with position and 
time; the second part, which arises from discrete encoun- 
ters with individual stars, fluctuates much more rapidly. 
The smooth force itself consists of the restoring force 
arising from the potential of the aggregate distribution 
of stars, and the dissipative force of dynamical friction 
which causes the object to decelerate as it moves through 
the stellar background . This is similar to the Langevin 
model of Brownian motion, which describes the irregu- 
lar motions of dust grains immersed in a gas of light 
molecules: the Brownian particle experiences a frictional 
force proportional to its velocity, and a random, rapidly 
fluctuating force owing to the large rate of collisions it 
suffers with the gas molecules in its neighborhood. 

In a stellar system, the analog of the gas molecules is 
the stars, and the Brownian particle corresponds to any 
object which is much more massive than the stars and 
moves much more slowly than they do. An astrophys- 
ically relevant example is a massive black hole at the 
center of a dense stellar system such as a globular cluster 
or the nucleus of a galaxy. We would like to extend the 
Langevin method of analysis to such a problem. 

Consider, therefore, a black hole of mass m in a cluster 
of stars of total mass M which we take to be described by 
the spherically symmetric density and potential profiles 
p(r) and <&(r), respectively, where r is the radial posi- 
tion vector from the center of the stellar system, which is 
taken as the origin. The phase space distribution func- 
tion, which in general will depend both on position r, 
and velocity v, is defined such that /(r, v) d 3 r d 3 v is the 
mass in stars in the phase space volume d 3 rd 3 v. We 
assume that in a spherically symmetric stellar model, / 



is a function only of the energy per unit mass E, where 
E= \v 2 + <&{r). 

A black hole in this stellar system is subject to three 
forces. The first is the restoring force of the stellar poten- 
tial, F = -mV$(r), where r is the position vector of the 
black hole. For a stellar system with a non-cuspy core, 
this reduces to the form F = —kr, where k is independent 
of r for small values of r. 

The second force on the black hole is the dissipative 
force of dynamical friction which causes it to decelerate 
as it moves through the stellar background. We use for 
this the well-known Chandrasekhar formula ( ||) F = 
— pV, where 



(3 = 167r 2 lnA G 2 m(m + m*) 



So /( r : u)u 2 du 



(1) 



In the above, v is the velocity of the black hole, m* 
is the mass of each star (in the following, we take all 
stars to have equal masses, for simplicity), and InA is the 
"Coulomb logarithm" ||. Since the black hole remains 
close to the origin and moves very slowly compared with 
the stars, we may replace f(r, u) in the integral by /(r, 0), 
and evaluate it at r — ► to obtain |H|: 



(3 -> (167r 2 /3)lnA G 2 m{m + m*)/(0, 0). 



(2) 



The third force on the black hole fluctuates on a time- 
scale which is extremely short compared to the above 
two forces [0, and arises from random discrete encoun- 
ters between the black hole and the stars. Denoting this 
stochastic force as F(t), we can write the equation of 
motion of the black hole as 



mr(t) + pr(i) + kr(t) = F(t), 



(3) 



which is the equation of motion of a harmonically bound 
Brownian particle. 
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The spatial components of this linear vector equation 
are separable into equivalent components: 



mx(t) + f3x(t) + kx{t) = F x (t). 



(4) 



The stochastic force is defined only statistically. Since it 
is random and rapidly varying, we expect it to be zero on 
average and uncorrelated with itself at different times: 

(F B (t))=0, {F x {t l )F x {t 2 ))=C8{t 1 -t 2 ) : (5) 

where 8 is the Dirac delta function and the angular brack- 
ets denote an average over an ensemble of "similarly pre- 
pared" systems of stars in each of which the black hole 
has the same initial position and velocity. We will later 
show how to determine the value of C. 

We have shown in detail in Ref. M how Eqs. (0) and 
(H) can be combined with the Fokker-Planck equation 
to derive the probability distributions of the black hole's 
position and velocity in the stationary state (i.e., when 
initial transients - which decay exponentially fast - die 
out and these distributions become time-independent); 
these turn out to be independent Gaussian distributions: 



W{x) = V^t/ttC luq m exp[— (2j/C)Lulm 2 x 2 ~ l 



W{v x ) = a/27/ttCto exp[-(2 7 /C)m 2 ?;2] 



(0) 



(7) 



where 7 = [3 /2m and uiq — y/k/m. Note that v x = x. 

It was shown in || (their Eq. [8-66]) that the rate of 
change of kinetic energy of the black hole is given by 



d(mv 2 /2)/dt = 167r 2 lnAG 2 mm* 



/(r, u)udu 



m/{m i ,v) / f(r 1 u)u 2 du 



(8) 



The first term describes the "heating" of the black hole 
by fluctuations in the stellar system, and the second term 
describes the "cooling" or dissipative effect of dynamical 
friction. Averaged over the ensemble in the stationary 
state, the two terms above should sum to zero. Again, 
setting f(r, u) in integrals of the form J Q ... f(r, u)du to 
f(r, 0), we can derive 

(« 2 ) = (3m t /m)| f(r,u)udu^/f(r,0)=3(v 2 x ), (9) 

since the three velocity components are equivalent and 
independent of each other. From Eqs. @ and (0), 



(x 2 ) = C/4jm 2 iu%, (v 2 x ) = C/4-fm 2 



(10) 



Combining Eqs. (g) and @, we obtain the following 
expressions for (x 2 ) and C: 



(x 2 ) = (m*/mo;g) ^ f(r,u)udu/f(r,0), (11) 



C = 47mm* J f(r, u)udu I f(r, 0). (12) 

We wish to have a measure for how far the black hole 
is from equipartition of kinetic energy between itself and 
the stars surrounding it in the core of the stellar system. 
Since the square of the stellar velocity dispersion is 



f{r, u) Attu 2 du If f(r, u) Ami 2 du, (13) 



we obtain by using Eq. (g): 



where 



n 



(v ) / v 2 — t) m+/r 



3 JcT /( r > u)udu / °° f(r, u)u 2 du 
/(r,0) / °°/(r, U )u 4 d M 



(14) 



(15) 



If r\ is evaluated in the limit r — > 0, it measures the 
deviation of the black hole from equipartition of kinetic 
energy between it and the stars in the core surrounding 
it. When 77 = 1, there is exact equipartition. We can 
take various stellar models and examine what the black 
hole's stationary state dynamics would be in each case. 

Maxwellian distribution. Suppose the phase space 
distribution of stars to be given by a Maxwellian dis- 
tribution with velocity dispersion a: 



f(r,u) oc exp(-w 2 /2a 2 ). 



(16) 



A simple calculation then gives rj — 1. This is not surpris- 
ing: the analogy in this case is with a Brownian particle 
immersed in a gas of molecules with a well-defined tem- 
perature H ; in the steady state, such a Brownian particle 
will be in exact equipartition with the molecules. 

The King model. The simplest model for which the 
distribution function is a pure Maxwellian is the "isother- 
mal sphere" , which is not a very physical description of 
actual stellar systems since it has infinite total mass. A 
commonly used alternative that resembles the isothermal 
sphere at small radii (where stars have large absolute val- 
ues of energy per unit mass) is the King model (see ||, 
1): 

f(E)= ( Pl (2na 2 )-i(e- E ^ 2 -I) if E < {yf) 
K ' \ if E > 0, v ' 

where p\ is independent of E. In Eq. (|l5|), the upper 
velocity limit of the integrals is now \j2a 2 ^ instead of 
infinity, where £ = — <I>(r)/er 2 . A calculation reveals that 



(l-i^i)(^ErfU/?)-£ 



1/2 



§e 3/2 ) 



Erf (v?) - e /2 - ie /2 - & 5/2 



(18) 



where Erf [z) is the error function defined as Erf (z) 
2/y / 7t L e~ l dt. For large values of £, 
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King model 



FIG. 1. A plot of 77 as a function of £ = — <3>(r)/cr 2 for a 
King model. The Brownian particle at the center of a King 
stellar model approaches exact equipartition as £ becomes 
large. 
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1 + (8£ 3/2 /i5\/^ - l)£ e ~« + • • • 0(e~ 2 «)ior£ > 1. 



Fig. 1 shows the relation between 77 and £. 77 starts at a 
value of approximately 1.75 and falls asymptotically to 1 
as £ rises. The King model goes over into the isothermal 
sphere as £ — > 00. The (mild) deviation from equiparti- 
tion is because of the deviation of the distribution func- 
tion from the Maxwellian form and, in particular, the 
existence of an upper velocity limit above which the dis- 
tribution function vanishes. Actual stellar systems are 
well described by King models with the central potential 
chosen such that — <£>(0)/<7 2 lies between roughly 3 and 
10 ( @, [§, @, §). As Fig. 1 shows, black holes at the 
centers of such systems will be very close to equipartition 
with the stars. 

Polytropic models. Consider next stellar systems de- 
scribed by a polytropic distribution function which takes 
the simple form 



/(£)« 



(-E)' 




if E < 
if E > 0. 



Simple calculations give 



r) = (n + 5/2)/(n+ 1) for n > — 1. 



(19) 



(20) 



Note that 7/ for a polytropic model is independent of po- 
sition r, unlike in the case of the King model. For large 
values of n, 77 is close to 1, but for small values of n, there 
can be significant deviations from equipartition. 

The potential profile of a polytropic system is a so- 
lution of the Lane-Emden equation 0, which can be 
solved in terms of elementary functions in only two cases: 
71 = 7/2 (this case was considered in 0]) and n = —1/2. 
We examine these two cases in more detail. 

Polytropic model with n = 7/2. When n = 7/2, we 
get the following functional forms for the density and 
potential profiles of the stellar system: 



p(r) 



Mia 2 
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47r (r 2 + a 2 ) 5 / 2 



$(r) = 



-GM 



where G is the gravitational constant, M is the total mass 
of the stellar system in stars, and a is a length parameter. 
This is the well-known Plummer model j|] which provides 
a good fit to some globular clusters. In this case, 77 = 4/3, 
which is mildly different from equipartition. luq = k/m = 
GMjo? for small r, and 



C = (8GM/9a)7mm*, 



(21) 



(a; 2 ) = 2a 2 TO*/9m, (tj 2 ) = 2GMm ir /9am. (22) 

Polytropic model with n = —1/2. When n — —1/2, 
we get the following functional forms for the density and 
potential profiles of the stellar system: 

/ \ _ f M/(Ait 2 a 2 ) x sin(r/a)/r if r < ira 
P[T) ~\0 if r > ira, 



$(r) 



—GMwa.(r/a)/'Kr if r < ira 
-GM/r + GM/ira if r > ira. 



In this case, 77 = 4, which is very different from equipar- 
tition. to 2 , = k/m, = GM/3ira 3 for small r, and 



C = (8GAf/7ra)7mm Jr , 



(23) 



(a; 2 ) = 6a 2 m*/m, (tj 2 ) = 2GMm i ./wam. (24) 



(r 2 + a 2 ) 1 / 2 ' 



The expressions for (x 2 ) in Eqs. (E2|) and ( p4[ ) agree 
with the result of Bahcall & Wolf |1C|] to within factors 
of order unity. 

Numerical simulations. In order to confirm our re- 
sults, we have performed N-body simulations, using suit- 
ably modified versions of the fourth-order integrators of 
Aarseth jLl| developed by Quinlan jl2| , for the two poly- 
tropic stellar systems described above. The techniques 
used are detailed in jl| and |l2j] . In both cases, we use 
units such that G = M = 1, and Coulomb forces are 
softened by e = 5 x 10~ 3 in order to prevent numerical 
divergences. Note that the polytropic initial conditions 
for the stellar system are not steady state solutions in the 
presence of the black hole; the initial transient response 
was therefore numerically settled. 

For the case of the Plummer model, we have taken the 
mass of the black hole to be m = 0.01, and the length 
scale of the potential to be a = 37r/16. The stellar system 
is made up of 100,000 stars, each of mass to* = 1 x 10~ 5 . 
The model was integrated for 600 time units. The left 
half of Fig. 2 shows the binned distributions of one com- 
ponent of the black hole's position and velocity in the 
stationary state. The solid lines show the corresponding 
bin values as predicted by Eqs. (|^) and (Q), under the as- 
sumption that in the stationary state, ensemble averages 
will be equal to time averages over long time spans. The 
agreement with the predictions of the model is evidently 
good (these results were described in B). 
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FIG. 2. The panels on the left show the distributions of 
x and v x from simulation (bins) and theory (line) for a black 
hole of mass m = 0.01 in a Plummer (n = 7/2 polytrope) 
model made up of 100,000 stars of total mass M = 1. The 
panels on the right show the corresponding results for a black 
hole of mass m = 0.03 in a n — —1/2 polytrope model of total 
mass M = 1 and 20,000 stars. 



non-singular at small radii; in realistic systems, however, 
the massive black hole would induce a density cusp con- 
sisting of stars bound tightly to it. Since it carries the 
cusp with it as it moves around, it is as if a black hole of 
a somewhat larger effective mass were moving in a back- 
ground of unbound stars whose density profile is fiat near 
the center (and unaffected by the black hole away from 
the center provided its mass is much smaller than the to- 
tal mass of the stellar system). Since the restoring force 
and dynamical friction are provided mainly by the un- 
bound stars, we would expect our results to apply to real 
systems as long as the mass of stars in the cusp is much 
smaller than the mass of the black hole; our simulations 
show that this is indeed so for the cases considered here. 
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The right half of Fig. 2 shows the corresponding 
plots for the n = —1/2 polytrope. We have here taken 
m = 3 x 10 ~ 2 , and a — 3/4ir. The stellar system is made 
up of 20,000 stars, each of mass = 5 x 10~ 5 . The 
integration time was 900 units. Agreement with the ana- 
lytical results is again good. We note that the agreement 
seen in Fig. 2 would be spoiled if the appropriate value 
for r\ were different from that given by our model. 

Summary. We conclude that the model outlined in 
this Letter provides a good statistical description of the 
dynamics of a Brownian particle - such as a massive black 
hole - in a stellar system. If the black hole is inserted at 
the center of such a system, the transient effects of its ini- 
tial position and velocity decay exponentially fast, and it 
settles into a stationary state in which the distributions 
of its position and velocity become time-independent. 
These distributions are also independent of each other 
and are Gaussian. If the stellar distribution function 
is Maxwellian, there is in the stationary state precise 
equipartition of kinetic energy between the black hole 
and the stars. The greater the deviation of the distribu- 
tion function from a Maxwellian form, the greater is the 
deviation from strict equipartition. We have shown sev- 
eral examples of such cases, in one of which the black hole 
is very far from equipartition. The deviation originates 
from the existence of an upper velocity limit for stars be- 
yond which the distribution function is zero. Note that 
the stationary state above is not a true equilibrium, since 
at very long time-scales (approximately 20 times the re- 
laxation time), the stellar system will undergo core col- 
lapse ||. 

The above analytical model applies strictly only to stel- 
lar models in which the density and potential profiles are 
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